In order of ascending contrast, we have found that the extended source fringes had the smallest contrast, followed by the point source fringes, followed then by the internal calibration source fringes. We do not take the etalon data into account here (which obviously have the largest contrast) given that an external etalon is used to produce the fringes. We are mainly interested by the fringes/interference caused by the MIRI detectors as such. As of yet it is not clear how the fringe contrast relates to the size of the source in the pupil plane. In fact if there exists such a trend, it is rather strange that the MTS extended source and the internal calibration (also extended) source do not look more similar in contrast. We saw however in Notebook 10 that the internal calibration source spectrum has spatial structure to it. This means that the source is resolved, or at least that it is less blurred out compared to the MTS extended source. Applying this to the fringes, these have more (or less) of an angle-averaging extent to them, thus they are less washed out for the internal calibration source compared to the MTS extended source. Let's see how this can be quantified.
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.''')
# import modules
import funcs
import mrsobs
import numpy as np
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
plt.style.use('presentation')
%matplotlib notebook
import warnings
warnings.simplefilter('ignore')
We load the images for one band of the MRS for different kinds of sources, including:
Additionally the pixel-to-wavelength calibration map and the pixel-to-along-slice position map are imported.
# Define paths to data
workDir = '/Users/ioannisa/Desktop/python/miri_devel/'
cdpDir = workDir+'cdp_data/'
d2cMapDir = workDir+'distortionMaps/'
lvl2path = workDir+'FM_data/LVL2/'
# Get data
band = '1A'
ext_source_sci,ext_source_bkg = mrsobs.FM_MTS_BB_extended_source(lvl2path,band,bb_temp='800K')
intcal_source_sci = mrsobs.MIRI_internal_calibration_source(lvl2path,band,campaign='FM')
point_source_sci_p1,point_source_bkg_p1 = mrsobs.FM_MTS_800K_BB_point_source_raster(lvl2path,position='middle',pointing='P3')
# Get wavelength calibration pixel map
d2cMaps = funcs.load_obj('d2cMaps_band{}'.format(band),path=d2cMapDir)
lambdaMap = d2cMaps['lambdaMap']
sliceMap = d2cMaps['sliceMap']
nslices = d2cMaps['nslices']
# Get spectral resolution table
specres_table = funcs.get_cdps(cdpDir=cdpDir,band=band)[4]
We subtract background exposures where available.
# perform transform
ext_source_bkgsubtr = ext_source_sci-ext_source_bkg
point_source_p1_bkgsubtr = point_source_sci_p1-point_source_bkg_p1
We perform an even-odd row signal correction to the data (caused by the read-out pattern of MIRI detector pixel rows).
ext_source_oddevencorr = funcs.OddEvenRowSignalCorrection(ext_source_bkgsubtr)
intcal_source_oddevencorr = funcs.OddEvenRowSignalCorrection(intcal_source_sci)
point_source_p1_oddevencorr = funcs.OddEvenRowSignalCorrection(point_source_p1_bkgsubtr)
One thing we do notice however is that the trend of increasing fringe amplitude with longer wavelengths is shared in all spectra. If so, could it be that the "only thing" differentiating the MTS extended source fringe transmission from the internal calibration source fringe (for instance) be a geometric or optical angle-averaging effect? Let us start where we left off in Notebook 10.
# pixel trace
ypos_p,xpos_p = funcs.detpixel_trace_compactsource(point_source_p1_oddevencorr,band,d2cMaps)
# signal normalization based on fringe profile
sel = (lambdaMap[ypos_p,xpos_p]!=0)
intcal_source_norm = funcs.norm_fringe(intcal_source_sci[ypos_p,xpos_p][sel],thres=0.,min_dist=6,k=3,ext=3)
ext_source_norm = funcs.norm_fringe(ext_source_oddevencorr[ypos_p,xpos_p][sel],thres=0.,min_dist=6,k=3,ext=3)
fig,axs = plt.subplots(3,1,figsize=(18,13))
axs[0].set_title('Internal Calibration Source spectrum')
axs[0].plot(lambdaMap[ypos_p,xpos_p][sel],intcal_source_norm[0])
axs[0].plot(lambdaMap[ypos_p,xpos_p][sel][intcal_source_norm[1]],intcal_source_norm[0][intcal_source_norm[1]],'ro')
axs[0].plot(lambdaMap[ypos_p,xpos_p][sel],intcal_source_norm[2])
axs[0].set_ylim(350,830)
axs[0].set_ylabel('Signal [DN/sec]')
for plot in range(1,3):
axs[plot].plot(lambdaMap[ypos_p,xpos_p][sel],intcal_source_norm[0]/intcal_source_norm[2],label='Internal Calibration Source')
axs[plot].plot(lambdaMap[ypos_p,xpos_p][sel],ext_source_norm[0]/ext_source_norm[2],label='MTS Extended Source')
axs[plot].hlines(1,4.88,5.77,linestyle='dashed')
axs[plot].set_ylim(0.7,1.05)
axs[plot].set_ylabel('Normalized signal')
axs[1].legend(loc='upper right',bbox_to_anchor=(0.65,1.08),fontsize=10)
axs[2].hlines([0.795,0.912],4.88,5.77,'r',alpha=0.4,linestyle='dashed')
axs[2].set_xlim(5.14,5.31)
axs[2].set_title('<Zoomed view>')
for plot in range(3): axs[plot].set_xlabel('Wavelength [micron]')
plt.tight_layout()
MTS extended source fringe scanning -->
# Compute zeroth order optical thickness of detector from fringe peak separation
ext_source_fringepeaks_wavelength = lambdaMap[ypos_p,xpos_p][ext_source_norm[1]] # microns
ext_source_fringepeaks_wavenumber = np.flipud(10000./ext_source_fringepeaks_wavelength) # cm-1
ext_source_mean_fringepeak_separation = np.mean(np.diff(ext_source_fringepeaks_wavenumber)[1:-1]) # omit first and last data point
ext_sel = np.nonzero(ext_source_norm[0])
ext_wavenumbers = np.flipud(10000./lambdaMap[ypos_p,xpos_p][ext_sel])
ext_norm_signal = np.flipud((ext_source_norm[0]/ext_source_norm[2])[ext_sel])
# define scanning window (arbitrarily set to 2*period)
N=1
while ext_wavenumbers[N]-ext_wavenumbers[0] < 2*ext_source_mean_fringepeak_separation: N+=1
# Scan fringes
ext_R,ext_D,ext_sigma_R,ext_sigma_D = [np.full(len(ext_wavenumbers),np.nan) for j in range(4)]
displacements = np.arange(-2,2,0.1) # um
len_numerics = []
for displacement in displacements:
R0,D0 = 0.10,1/(2*ext_source_mean_fringepeak_separation) + displacement/10000.
for i in range(len(ext_wavenumbers)-N):
wvnr_segment = ext_wavenumbers[i: N+i].copy()
array_segment = ext_norm_signal[i: N+i].copy()
# least-squares fitting
try:
popt,pcov = curve_fit(funcs.FPfunc_noPhaseShift,wvnr_segment,array_segment,p0=(R0,D0),absolute_sigma=True,method='lm')
ext_D[i: N+i] = popt[1]
except RuntimeError:
continue
len_numerics.append(len(np.where(np.abs(np.diff(ext_D))*10000.>1)[0]))
len_numerics = np.array(len_numerics)
R0,D0 = 0.10,1/(2*ext_source_mean_fringepeak_separation) + displacements[np.argmin(len_numerics)]/10000.
for i in range(len(ext_wavenumbers)-N):
wvnr_segment = ext_wavenumbers[i: N+i].copy()
array_segment = ext_norm_signal[i: N+i].copy()
# least-squares fitting
try:
popt,pcov = curve_fit(funcs.FPfunc_noPhaseShift,wvnr_segment,array_segment,p0=(R0,D0),absolute_sigma=True,method='lm')
TR = funcs.FPfunc_noPhaseShift(wvnr_segment,*popt)
ext_R[i: N+i] = popt[0]
ext_D[i: N+i] = popt[1]
ext_sigma_R[i: N+i] = np.sqrt(np.diag(pcov))[0]
ext_sigma_D[i: N+i] = np.sqrt(np.diag(pcov))[1]
except RuntimeError:
continue
# take care about numerics!
if len_numerics[np.argmin(len_numerics)] != 0:
ext_R,ext_D,numerics = funcs.cleanRD(ext_R,ext_D)
#> transmittance function
ext_TR = funcs.FPfunc_noPhaseShift(ext_wavenumbers,ext_R,ext_D)
# let's look at the results
fig,axs = plt.subplots(2,1,figsize=(18,8))
axs[0].plot(ext_wavenumbers,ext_R*100.)
axs[0].set_xlim(ext_wavenumbers[0],ext_wavenumbers[-1])
axs[0].set_ylim(1.2,4)
axs[0].set_ylabel('Reflectance [%]')
axs[1].plot(ext_wavenumbers,ext_D*10000.)
if len_numerics[np.argmin(len_numerics)] != 0:
axs[1].plot(ext_wavenumbers[numerics],ext_D[numerics]*10000,'ro',markersize=1,label='numerics')
axs[1].set_xlim(ext_wavenumbers[0],ext_wavenumbers[-1])
axs[1].set_ylabel('Optical thickness [micron]')
axs[1].set_xlabel('Wavenumber [cm-1]')
axs[1].legend()
plt.suptitle('Window size = {}'.format(N),fontsize=20)
plt.tight_layout(rect=[0, 0.03, 1, 0.96])
# let's apply the optimized transmittance function as our "fringe correction"
fig,axs = plt.subplots(3,1,figsize=(18,13))
axs[0].plot(ext_wavenumbers,ext_norm_signal,label='fringe transmission (in wavenumber space!)')
axs[0].plot(ext_wavenumbers,ext_TR,label='optimized transmittance function')
axs[0].hlines(1.05,ext_wavenumbers[0],ext_wavenumbers[-1],alpha=0.4,linestyle='dashed')
axs[0].plot(ext_wavenumbers,ext_norm_signal/ext_TR +0.05,label='residuals (offset by +0.05units)')
axs[0].hlines([1.04,1.06],ext_wavenumbers[0],ext_wavenumbers[-1],'r',alpha=0.4,linestyle='dashed',label='1% margin')
axs[0].set_xlim(ext_wavenumbers[0],ext_wavenumbers[-1])
axs[0].set_ylim(0.85,1.2)
axs[0].set_ylabel('FP Transmission')
axs[0].legend(loc='upper center')
axs[1].plot(ext_wavenumbers,np.flipud(ext_source_norm[0][ext_sel]),label='original spectrum')
axs[1].plot(ext_wavenumbers,np.flipud(ext_source_norm[2][ext_sel]),'g',label='spectral baseline')
axs[1].plot(ext_wavenumbers,np.flipud(ext_source_norm[2][ext_sel])*ext_TR,label='spectral baseline x optimized transmittance function')
axs[1].set_ylim(20,35)
axs[1].set_xlim(ext_wavenumbers[0],ext_wavenumbers[-1])
axs[1].set_ylabel('Signal [DN/s]')
axs[2].plot(ext_wavenumbers,np.flipud(ext_source_norm[0][ext_sel]),label='original spectrum')
axs[2].plot(ext_wavenumbers,np.flipud(ext_source_norm[0][ext_sel])/ext_TR,'r',label='fringe correction residuals')
axs[2].set_xlim(ext_wavenumbers[0],ext_wavenumbers[-1])
axs[2].set_ylim(20,35)
axs[2].set_ylabel('Signal [DN/s]')
for plot in range(1,3):
axs[plot].set_xlabel('Wavenumber [cm-1]')
axs[plot].legend(loc='upper left')
plt.tight_layout()
MIRI internal calibration source fringe scanning -->
# Compute zeroth order optical thickness of detector from fringe peak separation
intcal_source_fringepeaks_wavelength = lambdaMap[ypos_p,xpos_p][intcal_source_norm[1]] # microns
intcal_source_fringepeaks_wavenumber = np.flipud(10000./intcal_source_fringepeaks_wavelength) # cm-1
intcal_source_mean_fringepeak_separation = np.mean(np.diff(intcal_source_fringepeaks_wavenumber)[1:-1]) # omit first and last data point
intcal_sel = np.nonzero(intcal_source_norm[0])
intcal_wavenumbers = np.flipud(10000./lambdaMap[ypos_p,xpos_p][intcal_sel])
intcal_norm_signal = np.flipud((intcal_source_norm[0]/intcal_source_norm[2])[intcal_sel])
# define scanning window (arbitrarily set to 2*period)
N=1
while intcal_wavenumbers[N]-intcal_wavenumbers[0] < 2*intcal_source_mean_fringepeak_separation: N+=1
# Scan fringes
intcal_R,intcal_D,intcal_sigma_R,intcal_sigma_D = [np.full(len(intcal_wavenumbers),np.nan) for j in range(4)]
displacements = np.arange(-2,2,0.1) # um
len_numerics = []
for displacement in displacements:
R0,D0 = 0.10,1/(2*intcal_source_mean_fringepeak_separation) + displacement/10000.
for i in range(len(intcal_wavenumbers)-N):
wvnr_segment = intcal_wavenumbers[i: N+i].copy()
array_segment = intcal_norm_signal[i: N+i].copy()
# least-squares fitting
try:
popt,pcov = curve_fit(funcs.FPfunc_noPhaseShift,wvnr_segment,array_segment,p0=(R0,D0),absolute_sigma=True,method='lm')
intcal_D[i: N+i] = popt[1]
except RuntimeError:
continue
len_numerics.append(len(np.where(np.abs(np.diff(intcal_D))*10000.>1)[0]))
len_numerics = np.array(len_numerics)
R0,D0 = 0.10,1/(2*intcal_source_mean_fringepeak_separation) + displacements[np.argmin(len_numerics)]/10000.
for i in range(len(intcal_wavenumbers)-N):
wvnr_segment = intcal_wavenumbers[i: N+i].copy()
array_segment = intcal_norm_signal[i: N+i].copy()
# least-squares fitting
try:
popt,pcov = curve_fit(funcs.FPfunc_noPhaseShift,wvnr_segment,array_segment,p0=(R0,D0),absolute_sigma=True,method='lm')
TR = funcs.FPfunc_noPhaseShift(wvnr_segment,*popt)
intcal_R[i: N+i] = popt[0]
intcal_D[i: N+i] = popt[1]
intcal_sigma_R[i: N+i] = np.sqrt(np.diag(pcov))[0]
intcal_sigma_D[i: N+i] = np.sqrt(np.diag(pcov))[1]
except RuntimeError:
continue
# take care about numerics!
if len_numerics[np.argmin(len_numerics)] != 0:
intcal_R,intcal_D,numerics = funcs.cleanRD(intcal_R,intcal_D)
#> transmittance function
intcal_TR = funcs.FPfunc_noPhaseShift(intcal_wavenumbers,intcal_R,intcal_D)
# let's look at the results
fig,axs = plt.subplots(2,1,figsize=(18,8))
axs[0].plot(intcal_wavenumbers,intcal_R*100.)
axs[0].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
# axs[0].set_ylim(1.2,4)
axs[0].set_ylabel('Reflectance [%]')
axs[1].plot(intcal_wavenumbers,intcal_D*10000.)
if len_numerics[np.argmin(len_numerics)] != 0:
axs[1].plot(intcal_wavenumbers[numerics],intcal_D[numerics]*10000,'ro',markersize=1,label='numerics')
axs[1].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
axs[1].set_ylabel('Optical thickness [micron]')
axs[1].set_xlabel('Wavenumber [cm-1]')
axs[1].legend()
plt.suptitle('Window size = {}'.format(N),fontsize=20)
plt.tight_layout(rect=[0, 0.03, 1, 0.96])
# let's apply the optimized transmittance function as our "fringe correction"
fig,axs = plt.subplots(3,1,figsize=(18,13))
axs[0].plot(intcal_wavenumbers,intcal_norm_signal,label='fringe transmission (in wavenumber space!)')
axs[0].plot(intcal_wavenumbers,intcal_TR,label='optimized transmittance function')
axs[0].hlines(1.05,intcal_wavenumbers[0],intcal_wavenumbers[-1],alpha=0.4,linestyle='dashed')
axs[0].plot(intcal_wavenumbers,intcal_norm_signal/intcal_TR +0.05,label='residuals (offset by +0.05units)')
axs[0].hlines([1.02,1.08],intcal_wavenumbers[0],intcal_wavenumbers[-1],'r',alpha=0.4,linestyle='dashed',label='3% margin')
axs[0].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
axs[0].set_ylim(0.65,1.4)
axs[0].set_xlabel('Wavenumber [cm-1]')
axs[0].set_ylabel('FP Transmission')
axs[0].legend(loc='upper center')
axs[1].plot(intcal_wavenumbers,np.flipud(intcal_source_norm[0][intcal_sel]),label='original spectrum')
axs[1].plot(intcal_wavenumbers,np.flipud(intcal_source_norm[2][intcal_sel]),'g',label='spectral baseline')
axs[1].plot(intcal_wavenumbers,np.flipud(intcal_source_norm[2][intcal_sel])*intcal_TR,label='spectral baseline x optimized transmittance function')
axs[1].set_ylim(300,830)
axs[1].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
axs[1].set_ylabel('Signal [DN/s]')
axs[2].plot(intcal_wavenumbers,np.flipud(intcal_source_norm[0][intcal_sel]),label='original spectrum')
axs[2].plot(intcal_wavenumbers,np.flipud(intcal_source_norm[0][intcal_sel])/intcal_TR,'r',label='fringe correction residuals')
axs[2].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
axs[2].set_ylim(300,830)
axs[2].set_ylabel('Signal [DN/s]')
for plot in range(1,3):
axs[plot].set_xlabel('Wavenumber [cm-1]')
axs[plot].legend(loc='upper left')
plt.tight_layout()
Let's compare the optical properties derived from the MTS extended and the MIRI internal calibration source fringes.
fig,axs = plt.subplots(2,1,figsize=(18,8))
axs[0].plot(intcal_wavenumbers,intcal_R*100.,label='MIRI internal calibration source')
axs[0].plot(ext_wavenumbers,ext_R*100.,label='MTS extended source')
axs[0].set_ylabel('Reflectance [%]')
axs[1].plot(intcal_wavenumbers,intcal_D*10000.,label='MIRI internal calibration source')
axs[1].plot(ext_wavenumbers,ext_D*10000.,label='MTS extended source')
axs[1].set_ylabel('Optical thickness [micron]')
for plot in range(2):
axs[plot].set_xlim(intcal_wavenumbers[0],intcal_wavenumbers[-1])
axs[plot].set_xlabel('Wavenumber [cm-1]')
axs[plot].legend()
plt.tight_layout(rect=[0, 0.03, 1, 0.96])
We know that the wavelength solution is consistent from analyzing etalon data. As such, the discrepancy in the optical thickness must be caused by a geometric effect on the optics of the MIRI internal calibration source illumination. Based on that conjecture we use the optical thickness determined from the MTS extended source fringes as the optical thickness of the MIRI internal calibration source fringe transmission.
intcal_detector_R,intcal_detector_D = np.full(1024,np.nan),np.full(1024,np.nan)
ext_detector_R,ext_detector_D = np.full(1024,np.nan),np.full(1024,np.nan)
intcal_detector_R[intcal_sel],intcal_detector_D[intcal_sel] = np.flipud(intcal_R),np.flipud(intcal_D)
ext_detector_R[ext_sel],ext_detector_D[ext_sel] = np.flipud(ext_R),np.flipud(ext_D)
#> transmittance function
wavenumbers = np.flipud(10000./lambdaMap[ypos_p,xpos_p])
intcal_detector_TR = funcs.FPfunc_noPhaseShift(wavenumbers,np.flipud(intcal_detector_R),np.flipud(ext_detector_D) )
ext_detector_TR = funcs.FPfunc_noPhaseShift(wavenumbers,np.flipud(ext_detector_R),np.flipud(ext_detector_D) )
fig,axs = plt.subplots(2,1,figsize=(18,8))
axs[0].set_title('MIRI Internal Calibration Source')
for plot in range(2):
axs[plot].plot(intcal_wavenumbers,intcal_TR,label='optical thickness of internal calibration source')
axs[plot].plot(wavenumbers,intcal_detector_TR,label='optical thickness of extended source')
axs[plot].set_xlabel('Wavenumber [cm-1]')
axs[plot].set_ylabel('Transmittance function')
axs[1].set_xlim(1732.7,1816.5)
axs[1].set_title('<Zoomed view>')
axs[0].legend(loc='lower right')
plt.tight_layout()
The contrast of the internal calibration source fringes is much higher than that of the extended source fringes. Contrast can be reduced when averaging the fringe transmission over a range of phases. Let's see how this impacts the derived transmittance fuction of the internal calibration source fringes.
To apply the phase-averaging we:
After doing so we normalize the transmittance functions to the fringe peaks and compare the results.
fig,axs = plt.subplots(2,1,figsize=(18,8))
for plot in range(2):
axs[plot].plot(wavenumbers,intcal_detector_TR,label=r'$\beta = 0\pi$')
for beta in [np.pi/5,np.pi/4,np.pi/3,np.pi/2.]:
intcal_detector_TR_phase_avg = 0
for rand in np.random.uniform(-1.,1.,5000):
intcal_detector_TR_phase_avg += funcs.FPfunc(wavenumbers,np.flipud(intcal_detector_R),np.flipud(ext_detector_D),np.pi + beta*rand)
intcal_detector_TR_phase_avg /= 5000.
axs[0].plot(wavenumbers,intcal_detector_TR_phase_avg,label=r'$\beta = {} \pi$'.format(round(beta/np.pi,2)))
# normalize fringes
intcal_detector_TR_phase_avg_norm = funcs.norm_fringe(intcal_detector_TR_phase_avg,thres=0.,min_dist=2,k=3,ext=3)
axs[1].plot(wavenumbers,intcal_detector_TR_phase_avg_norm[0]/intcal_detector_TR_phase_avg_norm[2],label=r'$\beta = {} \pi$'.format(round(beta/np.pi,2)))
for plot in range(2):
axs[plot].hlines(1,1730,2050,linestyle='dashed')
axs[plot].set_ylim(0.70)
axs[plot].legend(loc='lower right')
axs[plot].set_xlabel('Wavenumber')
axs[plot].set_ylabel('Transmittance function')
plt.tight_layout()
It is not clear if one, two, three, or all four reasons above have an appreciable effect in the case of MIRI and the ground test observations. But let's compare the phase-averaged internal calibration source transmittance function to that of the MTS extended source. A $\chi^2$ analysis minimization is carried out to determine the value of $\beta$ that minimizes the difference between the two.
# determine partial coherence beta factor by minimizing difference between angle-averaged FP model and extended source data
betas = np.linspace(0,np.pi/2.,200)
# perform chi-square minimization
chi_square = []
for beta in betas:
intcal_detector_TR_phase_avg = 0
for rand in np.random.uniform(-1.,1.,1000):
intcal_detector_TR_phase_avg += funcs.FPfunc(wavenumbers,np.flipud(intcal_detector_R),np.flipud(ext_detector_D),np.pi + beta*rand)
intcal_detector_TR_phase_avg /= 1000.
intcal_detector_TR_phase_avg_norm = funcs.norm_fringe(intcal_detector_TR_phase_avg,thres=0.,min_dist=2,k=3,ext=3)
chi_square.append(np.sum(ext_detector_TR[~np.isnan(ext_detector_TR)]/(intcal_detector_TR_phase_avg_norm[0]/intcal_detector_TR_phase_avg_norm[2])[~np.isnan(ext_detector_TR)] -1 )**2)
beta = betas[np.argmin(np.array(chi_square))]
print r'Minimum chi^2 = {}'.format(round(np.min(np.array(chi_square)),5) )
print r'chi^2 minimized for beta = {} rad'.format(round(beta,4)) # radians
print r'chi^2 minimized for beta = {} pi rad'.format(round(beta/np.pi,2)) # radians
print r'chi^2 minimized for beta = {} deg'.format(round(beta*180./np.pi,4)) # degrees
plt.figure(figsize=(18,4))
plt.plot(betas/np.pi,chi_square)
plt.plot(beta/np.pi,np.min(np.array(chi_square)),'ro')
plt.hlines(np.min(np.array(chi_square)),0,0.5,linestyle='dashed')
plt.vlines(beta/np.pi,-200,4650,linestyle='dashed')
plt.ylim(-200)
plt.xlabel(r'$\beta$ [$\pi$ rad]',fontsize=20)
plt.ylabel(r'$\chi^2$',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
intcal_detector_TR_phase_avg = 0
for rand in np.random.uniform(-1.,1.,5000):
intcal_detector_TR_phase_avg += funcs.FPfunc(wavenumbers,np.flipud(intcal_detector_R),np.flipud(ext_detector_D),np.pi + beta*rand)
intcal_detector_TR_phase_avg /= 5000.
# normalize fringes
intcal_detector_TR_phase_avg_norm = funcs.norm_fringe(intcal_detector_TR_phase_avg,thres=0.,min_dist=2,k=3,ext=3)
plt.figure(figsize=(18,4))
plt.plot(wavenumbers,(intcal_detector_TR_phase_avg_norm[0]/intcal_detector_TR_phase_avg_norm[2]),label=r'Phase-averaged int cal source transm function ($\beta$={}$\pi$ rad)'.format(round(beta/np.pi,2) ))
plt.plot(wavenumbers,ext_detector_TR,label='MTS extended source transmittance function')
plt.plot(wavenumbers,ext_detector_TR/(intcal_detector_TR_phase_avg_norm[0]/intcal_detector_TR_phase_avg_norm[2]),label='residuals')
plt.hlines([0.99,1.01],1730,2050,'r',linestyle='dashed',alpha=0.4,label='1% margin')
plt.ylim(0.85,1.03)
plt.xlabel('Wavenumber [cm-1]')
plt.ylabel('Transmittance function')
plt.legend(fontsize=10)
plt.tight_layout()
We find that phase-averaging the transmittance function of the MIRI internal calibration source results in a better match between the int.cal.source fringes and that of the MTS external calibration source. This phase-averaging process can be used as a tool to find systematics in the comparison between the int.cal.source fringes and that of the external sources.
Can we apply the collection of analysis techniques developed to study global variations of parameters over the detector plane? In the notebooks that follow we perform large-scale analyses of fringe data to understand how optical parameters vary: